home *** CD-ROM | disk | FTP | other *** search
/ Enter 2004 January / enter-2004-01.iso / files / maxima-5.9.0.exe / {app} / share / maxima / 5.9.0 / src / numerical / slatec / zbesy.lisp < prev    next >
Encoding:
Text File  |  2003-02-09  |  4.5 KB  |  122 lines
  1. ;;; Compiled by f2cl version 2.0 beta 2002-05-06
  2. ;;; 
  3. ;;; Options: ((:prune-labels nil) (:auto-save t) (:relaxed-array-decls t)
  4. ;;;           (:coerce-assigns :as-needed) (:array-type ':simple-array)
  5. ;;;           (:array-slicing nil) (:declare-common nil)
  6. ;;;           (:float-format double-float))
  7.  
  8. (in-package "SLATEC")
  9.  
  10.  
  11. (defun zbesy (zr zi fnu kode n cyr cyi nz cwrkr cwrki ierr)
  12.   (declare (type (simple-array double-float (*)) cwrki cwrkr cyi cyr)
  13.            (type f2cl-lib:integer4 ierr nz n kode)
  14.            (type double-float fnu zi zr))
  15.   (prog ((i 0) (k 0) (k1 0) (k2 0) (nz1 0) (nz2 0) (c1i 0.0) (c1r 0.0)
  16.          (c2i 0.0) (c2r 0.0) (elim 0.0) (exi 0.0) (exr 0.0) (ey 0.0) (hcii 0.0)
  17.          (sti 0.0) (str 0.0) (tay 0.0) (ascle 0.0) (rtol 0.0) (atol 0.0)
  18.          (aa 0.0) (bb 0.0) (tol 0.0) (r1m5 0.0))
  19.     (declare
  20.      (type double-float r1m5 tol bb aa atol rtol ascle tay str sti hcii ey exr
  21.       exi elim c2r c2i c1r c1i)
  22.      (type f2cl-lib:integer4 nz2 nz1 k2 k1 k i))
  23.     (setf ierr 0)
  24.     (setf nz 0)
  25.     (if (and (= zr 0.0) (= zi 0.0)) (setf ierr 1))
  26.     (if (< fnu 0.0) (setf ierr 1))
  27.     (if (or (< kode 1) (> kode 2)) (setf ierr 1))
  28.     (if (< n 1) (setf ierr 1))
  29.     (if (/= ierr 0) (go end_label))
  30.     (setf hcii 0.5)
  31.     (multiple-value-bind
  32.         (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9)
  33.         (zbesh zr zi fnu kode 1 n cyr cyi nz1 ierr)
  34.       (declare (ignore var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7))
  35.       (setf nz1 var-8)
  36.       (setf ierr var-9))
  37.     (if (and (/= ierr 0) (/= ierr 3)) (go label170))
  38.     (multiple-value-bind
  39.         (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9)
  40.         (zbesh zr zi fnu kode 2 n cwrkr cwrki nz2 ierr)
  41.       (declare (ignore var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7))
  42.       (setf nz2 var-8)
  43.       (setf ierr var-9))
  44.     (if (and (/= ierr 0) (/= ierr 3)) (go label170))
  45.     (setf nz (min (the f2cl-lib:integer4 nz1) (the f2cl-lib:integer4 nz2)))
  46.     (if (= kode 2) (go label60))
  47.     (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
  48.                   ((> i n) nil)
  49.       (tagbody
  50.         (setf str
  51.                 (- (f2cl-lib:fref cwrkr (i) ((1 n)))
  52.                    (f2cl-lib:fref cyr (i) ((1 n)))))
  53.         (setf sti
  54.                 (- (f2cl-lib:fref cwrki (i) ((1 n)))
  55.                    (f2cl-lib:fref cyi (i) ((1 n)))))
  56.         (f2cl-lib:fset (f2cl-lib:fref cyr (i) ((1 n))) (* (- sti) hcii))
  57.         (f2cl-lib:fset (f2cl-lib:fref cyi (i) ((1 n))) (* str hcii))
  58.        label50))
  59.     (go end_label)
  60.    label60
  61.     (setf tol (max (f2cl-lib:d1mach 4) 1.0e-18))
  62.     (setf k1 (f2cl-lib:i1mach 15))
  63.     (setf k2 (f2cl-lib:i1mach 16))
  64.     (setf k (f2cl-lib:int (min (abs k1) (abs k2))))
  65.     (setf r1m5 (f2cl-lib:d1mach 5))
  66.     (setf elim (* 2.303 (- (* k r1m5) 3.0)))
  67.     (setf exr (cos zr))
  68.     (setf exi (sin zr))
  69.     (setf ey 0.0)
  70.     (setf tay (coerce (abs (+ zi zi)) 'double-float))
  71.     (if (< tay elim) (setf ey (exp (- tay))))
  72.     (if (< zi 0.0) (go label90))
  73.     (setf c1r (* exr ey))
  74.     (setf c1i (* exi ey))
  75.     (setf c2r exr)
  76.     (setf c2i (- exi))
  77.    label70
  78.     (setf nz 0)
  79.     (setf rtol (/ 1.0 tol))
  80.     (setf ascle (* (f2cl-lib:d1mach 1) rtol 1000.0))
  81.     (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
  82.                   ((> i n) nil)
  83.       (tagbody
  84.         (setf aa (f2cl-lib:fref cwrkr (i) ((1 n))))
  85.         (setf bb (f2cl-lib:fref cwrki (i) ((1 n))))
  86.         (setf atol 1.0)
  87.         (if (> (max (abs aa) (abs bb)) ascle) (go label75))
  88.         (setf aa (* aa rtol))
  89.         (setf bb (* bb rtol))
  90.         (setf atol tol)
  91.        label75
  92.         (setf str (* (- (* aa c2r) (* bb c2i)) atol))
  93.         (setf sti (* (+ (* aa c2i) (* bb c2r)) atol))
  94.         (setf aa (f2cl-lib:fref cyr (i) ((1 n))))
  95.         (setf bb (f2cl-lib:fref cyi (i) ((1 n))))
  96.         (setf atol 1.0)
  97.         (if (> (max (abs aa) (abs bb)) ascle) (go label85))
  98.         (setf aa (* aa rtol))
  99.         (setf bb (* bb rtol))
  100.         (setf atol tol)
  101.        label85
  102.         (setf str (- str (* (+ (* aa c1r) (* -1 bb c1i)) atol)))
  103.         (setf sti (- sti (* (+ (* aa c1i) (* bb c1r)) atol)))
  104.         (f2cl-lib:fset (f2cl-lib:fref cyr (i) ((1 n))) (* (- sti) hcii))
  105.         (f2cl-lib:fset (f2cl-lib:fref cyi (i) ((1 n))) (* str hcii))
  106.         (if (and (= str 0.0) (= sti 0.0) (= ey 0.0))
  107.             (setf nz (f2cl-lib:int-add nz 1)))
  108.        label80))
  109.     (go end_label)
  110.    label90
  111.     (setf c1r exr)
  112.     (setf c1i exi)
  113.     (setf c2r (* exr ey))
  114.     (setf c2i (* (- exi) ey))
  115.     (go label70)
  116.    label170
  117.     (setf nz 0)
  118.     (go end_label)
  119.    end_label
  120.     (return (values nil nil nil nil nil nil nil nz nil nil ierr))))
  121.  
  122.